%% MICRO-428, Week Eight: HDR Exercise Solution

% Build an High Dynamic Range (HDR) Image by stacking 4 different
% exposure images. Use HSL colour space for simple Lightness manipulation
% of the four images.

clear all;
close all;
clc;

figure;
% Load figures in decreasing brightness order.
imgs = {'Expo7.jpg','Expo5.jpg','Expo3.jpg','Expo1.jpg'};
%expTimes = [1/60, 1/120,1/250,1/500,1/1000,1/2000,1/4000]; % Exposure
%Times seven images
expTimes = [1/60,1/250,1/1000,1/4000];
montage(imgs) % Make a montage of the pics

% Get the images number:
N = length(imgs); % Count images

% Get the number of pixels of the pictures:
photo1 = imread('Expo7.jpg'); 
% Get image size
SizePhoto = []; 
% Get total number of pixels
PixelSize = []; 

%% Plot a Histogram for each Picture

for i =1:N
    img = imread(imgs{i});
    figure('Name',['Picture Exposure # ',num2str(i)])
    subplot(1,2,1);
    imshow(img)
    subplot(1,2,2);
    makePictureHistogram(img);
end

%% Test the pictures with matlab native HDR function

% Add makehdr function here. 
hdr =[];

% Plot the hdr image
figure('Name',' MATLAB HDR Image');
subplot(1,2,1);
rgb = tonemap(hdr);
imshow(rgb)
subplot(1,2,2);
makePictureHistogram(hdr);

% Save the image on the main folder
imwrite(rgb,'Img_MATLAB_HDR.jpg')

%% Set underexposed and overexposed values: analysis on R, G, B
 
% Set to allocate memory for the HDR photo, the properly exposed counter,
% and the positions of properly exposed, underexposed and overexposed 
% pixels.
hdrImageHsl = [];
properlyExposedCount = [];
pixelUnderExposed = [];
pixelOverExposed = [];


% For te weights, we could inn principle use the relative exposure values:
relW = [];
% However, in case we didnt know the exposure time or the ISO change, we
% could measure the exposure ratios. Let's prepare an empty array that we
% will fill later on with the measured values.
measW = zeros(1,N);
 
% Now, we have to follow the algorithm.

% Set the lower and higher limits relative to overexposition and
% underexposition
% Add a lower bound for noise
LowerLimit=[];
% Add an higer bound for saturation
UpperLimit =[];
% Estimate the forced maximum dynamic range of the images:
StartDynRange = UpperLimit/LowerLimit;

% Final preparation component we need before we jump in the algorithm 
% itself, is a measure of the intensity. Let's assume we do not know the
% sutter speeds, but we rather extrapolate them from the pictures.

% First, let's find the non-saturated pixels common to both the brightest  
% and the darkest images of the stack: they are going to be used as sample
% pixels to estimate the overall exposure value of the entire picture.
imgOE = rgb2hsl(double(imread(imgs{1})));
imgUE = rgb2hsl(double(imread(imgs{N})));
% Find the underexposed pixels. 
OEunderExposed = [];
% Find the overexposed pixels.
OEoverExposed = [];
% Find the properly exposed pixels.
ProperlyExposed = ~(OEunderExposed | OEoverExposed);

% Let's define the maximum intensity of our scene by measuring the
% intensity of the overexposed image.
intensity0 = [];

%% Main algorithm

for i = 1:N
    % Read the data of the image (RGB)
    img = double(imread(imgs{i}));
    % convert it from RGB format to HSL
    imgHsl = rgb2hsl(img);
    % Split H, S and L array for ease of manipulation
    imgH = imgHsl(:,:,1);
    imgS = imgHsl(:,:,2);
    imgL = imgHsl(:,:,3);
        
    % Display Dynamic Range
    DynRange = max(max(imgL))/max(min(min(imgL)),LowerLimit);
    disp(['Image ',num2str(i),': Dynamic Range = ',num2str(DynRange)]);
    
    % The following algorithm will act on the Luminosity value of the HSL
    % transform of the RGB file.

    % Find the underexposed pixels. Update the underexposed pixel matrix
    underExposed = imgL < LowerLimit;
 
    % Find the overexposed pixels. Update the overexposed pixel
    overExposed = imgL > UpperLimit;
    
    % Find the properly exposed pixels. Keep the count
    properlyExposed = ~(underExposed | overExposed);
    AllproperlyExposed = repmat(properlyExposed,1,1,3);
    properlyExposedCount(AllproperlyExposed) = properlyExposedCount(AllproperlyExposed) + 1;
    
    % First, let's measure the Bright and Dark pixels percentage of each
    % figure.
    DarkPerc = sum(sum(double(imgL <= LowerLimit)))/PixelSize;
    BrightPerc = sum(sum(double(imgL >= UpperLimit)))/PixelSize;
    disp(['Dark Pixels Percentage: ',num2str(DarkPerc*100,'%3.2f'),'%'])
    disp(['Bright Pixels Percentage: ',num2str(BrightPerc*100,'%3.2f'),'%'])
    disp(' ')
    
    % Let's measure the intensity of the picture
   intensity = sum(imgL(:,:));
    % Let's allocate the value of the relative intensity
    measW(i) = intensity0/intensity;
    
    % Set to zero all the underexposed and overexposed pixels.
    imgL(~properlyExposed) = 0;
    imgH(~properlyExposed) = 0;
    imgS(~properlyExposed) = 0;
    
    % Decide which weights to use (relW or measW):
    weights = relW;
    
    % Update the HDR image
    hdrImageHsl(:,:,1) = hdrImageHsl(:,:,1) + imgH;
    hdrImageHsl(:,:,2) = hdrImageHsl(:,:,2) +imgS;
    hdrImageHsl(:,:,3) = [];

    
    % Before going on, let's make a visual image for overexposed, 
    % underexposed and properly exposed pixels.
    imgExpo = zeros(SizePhoto(1),SizePhoto(2));
    imgExpo = imgExpo + 200.0 * ones(SizePhoto(1),SizePhoto(2 ));
    imgExpo = imgExpo - 200.0 * double(underExposed);
    imgExpo = imgExpo + 55 * double(overExposed);
    imwrite(uint8(imgExpo),['Img_Expo',num2str(i),'.jpg']);
    
end
%plot helper images
figure()
 imgs_Expo =flip( {'Img_Expo7.jpg','Img_Expo5.jpg','Img_Expo3.jpg','Img_Expo1.jpg'});
 montage(imgs_Expo)

%% Conversion back to RGB 

% Divide the pixels that we counted for more than one time by the number
% of times we counted them
hdrImageHsl = [];

% Rescale the luminosity L value back to the range [0,UpperLimit]
hdrImageHsl(:,:,3) = hdrImageHsl(:,:,3) / max(max(hdrImageHsl(:,:,3))) * UpperLimit;
% Set the overexposed pixels of the max exposure pic to L=1
% Add the *unproper* pixel values here.

UEunderExposed = imgUE(:,:,3) < LowerLimit;
UEoverExposed = imgUE(:,:,3) > UpperLimit;
hdrImageHsl = hdrImageHsl + double(UEoverExposed);
hdrImageHsl = hdrImageHsl + double(OEunderExposed);

% Measure Final Dynamic Range
FinalDynRange = max(max(hdrImageHsl(:,:,3)))/min(min(hdrImageHsl(:,:,3)));
disp(['HDR Image: Dynamic Range = ',num2str(FinalDynRange)]);

% Convert image back from HSL to RGB
hdrImage = hsl2rgb(hdrImageHsl);
    
figure('Name','Custom HDR Image');
subplot(1,2,1);
hdr = tonemap(hdrImage);
imshow(hdr)
subplot(1,2,2);
makePictureHistogram(hdrImage);

imwrite(hdr,'Img_HDR_HSL.jpg')
